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Abstract 

The intensity distribution of electromagnetic polar waves in a chain of near-resonant weakly- 
coupled scatterers is investigated theoretically and supported by a numerical analysis. Critical 
scaling behavior is discovered for part of the eigenvalue spectrum due to the disorder-induced An- 
derson transition. This localization transition (in a formally one-dimensional system) is attributed 
to the long-range dipole-dipole interaction, which decays inverse linearly with distance for polar- 
ization perpendicular to the chain. For polarization parallel to the chain, with inverse squared 
long range coupling, all eigenmodes are shown to be localized. A comparison with the results for 
Hermitian power-law banded random matrices and other intermediate models is presented. This 
comparison reveals the significance of non-Hermiticity of the model and the periodic modulation 
of the coupling. 
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I. INTRODUCTION 



Collective excitations of nanoparticle composites have shown promising applications for 
sensing, nonlinear spectroscopy, and photonic circuits. Among these applications, transport 
of electromagnetic signal along an assembly of metallic nanoparticles has been the subject 
of intensive research in recent years. It has shown promising applications in integrated 
photonics^ and sensing^. By the nature of their fabrication, disorder is inevitable in these 
artificial structures and therefore must be considered accordingly. 

In this article, we make a connection between the photonic transport in these novel 
physical structures and the Anderson localization transition: a fundamental phenomenon 
which emerges in various fields such as condensed matter physics^, cold gases in optical 
lattices^ and classical waves in random media^. We argue how the polar excitations in a 
chain of resonators can show critical scaling behavior. This criticality plays a major role in 
understanding the underlying phase transition phenomena. 

Anderson localization in electronic systems is extensively studied in the form of the 
Anderson tight-binding Hamiltonian with on-site disorder-. For this model all states are 
exponentially-localized in one and two dimensions. In three dimensional space there exists a 
metal-insulator transition, as suggested by the single parameter scaling theory. At the An- 
derson transition, wavefunctions show critical statistical behavior and their spatial structure 
are multifractaF'^'*. 

Levitov^*^ has studied the effect of long-range interaction, in the form of Dij oc |ri — rj|~'^, 
on the the localization of vibrational excitations of a disordered lattice. He showed that for 
fi > d, where d is the dimensionality of space, all states remain localized. For yU < ci all 
states escape localization due to a diverging number of resonances between spectrally apart 
energy levels. For fi = d, delocalization is weak and states are critical. 

Transition from localized to extended eigenstates has been shown by Mirlin et. al.— in 
the ensemble of power- law random banded matrices (PLRBM). In these Hermitian random 
matrices the off-diagonal elements decay as = {b/\i — for \i — j\ > b, where b is 

called the bandwidth. Similar to Levitov's model all eigenvectors are localized for /i > 1, the 
tight-binding limit, and are extended (metallic) for /i < 1, as for conventional Wigner-Dyson 
random matrices. At yU = 1, eigenstates show critical statistics for any width of the band. 

A localization transition is also shown for a non-Hermitian disordered tight-binding model 
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by Hatano and Nelsoiii^. Their model was motivated by its application to a special mapping 
of flux lines in certain superconductors to a bosonic system with a random potential. 

The physical system considered in this report, is fully described by a class of complex- 
symmetric coupling matrices. The eigenvectors of these matrices describe the excitation 
"modes" of a linear chain of point-like scatterers. These scatterers are driven close to 
resonance and are coupled to each other through long-range dipole-dipole interaction. 

We have studied this system in two cases of weak and strong coupling. By studying the 
scaling behavior of eigenstates in the weak coupling regime, we show that for transverse 
magnetic (TM) polarization parallel to the chain direction, all the states are localized. For 
the transverse electromagnetic (TEM) polarization in this regime, some of the states are 
critically extended and their scaling is described by a multifractal spectrum. We analyt- 
ically derive a perturbation expression for this multifractal spectrum in the limit of weak 
coupling corresponding to very strong disorder. In the strong coupling regime, we show 
some numerical evidence that the intensity distribution follows a mixed phase of localized 
and extended statistics. 

We have extensively compared the scaling behavior of this physical system with several 
hypothetic Hermitian and non-Hermitian matrix ensembles. This comparison proves the 
strong influence of the phase periodicity in the coupling terms. For example, the Levitov 
matrices with fi = 1 are no longer critical, but localizing, if the off-diagonal terms are 
non-random and have a periodic phase relation. On the other hand, the observed critical 
behavior of TEM polar eigenmodes disappear if the interaction phase factor is randomized. 

Our findings provide a clear and universal framework for excitation properties of an 
important building block in modern photonics. On a broader perspective our model has 
significant resemblance with the other important classes of Hamiltonians, which are used 
for describing several transport phenomena in mesoscopic systems. Since our model has 
an exact correspondence to a real physical system, it will pave the way for experimental 
investigation of several theoretical findings, which up to now were bound to the limitations 
of numerical simulation. 
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II. THE MODEL 

For describing the chain of resonators, we use the dipole approximation for each of the 
scatterers and the full dyadic on-shell Green function for their interaction. This model is 
previously used for describing collective plasmon excitations of metallic nanoparticles on 
a line or a plane for periodio^*^, aperiodic^ and disordered configurations^^. In partic- 
ular Markel and Sarychev have reported signatures of localization in a chain of point-like 
scatterersi^. 

The presence of an Anderson transition, its critical behavior, and the detailed statistics 
of localized or delocalized modes in such a system has not yet been studied. In the following 
we will argue and show analytically that a disorder mediated delocalization transition can 
happen for polarization perpendicular to the chain direction (TEM modes), while for po- 
larization parallel to the chain direction (TM modes), all eigenstates are localized in a long 
enough chain. We will present our results using a well-established statistical framework of 
probability density function (PDF) of eigenmode intensities and the scaling of generalized 
inverse participation ratios (GIPR). 

A. Dipole chain model 

We consider a linear array of L equally spaced polarizable isotropic particles with an 
interparticle distance of s. The size of the particles are considered small enough, relative 
to both s and the excitation wavelength A = 2Tic/uj, for a point-dipole approximation to 
be valid. With these considerations, TEM and TM modes are decoupled from each other. 
For the stationary response, oscillating with constant frequency w, the dipole moments of 
particles Pi = u.pj projected on each mode, are the solutions to the following homogeneous 
set of linear equations 



where is the polarizability of the ith particle and e*'^*-E™'^(a;j) is the incident electric field 
at its position projected on the specific Cartesian coordinate of the mode, u. The free space 
Green function should be replaced by the proper expressions for TEM(±) and TM(||) 




(1) 
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modes, which are given by 



^"^^ 47re Vc% ca;2 



(2) 




(3) 



Mij = 5ija^ + (1 - 5ij)g^{\xi - Xj\). 



(4) 



The exphcit frequency dependence of is dropped, since we consider only monochromatic 
excitations in this article. The matrix is a complex and symmetric matrix, the inverse 
of which gives the polarization response of the system to an arbitrary excitation; \p) = 



points. Since is non- Hermit ian, its eigenvalues are complex. However, the eigenvectors 
form a complete (bi-orthogonal) basis, unless the matrix is defective. Defective matrices 
form a subset of measure zero for a randomly generated ensemble. No defective matrix has 
been encountered in this research. The orthogonality condition is set by the quasi-scalar 
product of each two eigenvectors: 



where {ipn) is a right eigenvector of A^; Ai \ipn) = iV'n)- The eigenvectors are normalized 
to unity: {ipnl '4'n) = 1- The quasi-scalar product of an eigenvector with itself is a non- zero 
complex number for non-defective matrices. 

Under the stated assumptions, the polarization response to an incident field can be ob- 
tained from the decomposition 



A null eigenvalue points to a collective resonance of the system and the corresponding 
eigenvector is the most bound (guided) mode with the highest polarizability. 



Ai ^ \E^^^y In fact ^ is the so-called t-matrix^^ of the chain specified on the lattice 




(5) 




(6) 
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B. Resonant point scatterer 



A simple and yet general model for the dipolar polarizability of a point scatterer that 
conserves energy^ is given by 



where the last term in Eq. is the first non-vanishing radiative correction that fullfils the 
optical theorem. The quasistatic polarizability a° depends on the particle shape and its 
material properties. For a Lorentzian resonance around ur 



where V is the volume of the scatterer, 7 is the Ohmic damping factor and A is a constant 
that depends only on the geometry of the scatterer. For elastic scatterers a° is real-valued 
and diverges on resonance. 

C. Dimensionless formulation 

To study the properties of the coupling matrix (jl]) both theoretically and numerically, we 
rewrite it in terms of dimensionless quantities by dividing all the length dimensions by the 
interparticle distance s and multiplying the unit of polarizability by Airek^, where k = u/c. 
For the cases considered in this article, we also neglect the Ohmic damping of scatterers and 
hence the imaginary part on the diagonal of the matrix is given by the radiative damping 
term in Eq. ([7]). 

Based on definition (jl]) two distinct types of disorder can be considered for the system 
under investigation. Pure off-diagonal disorder is caused by the variation in the inter-particle 
spacing considering identical scatterers. The contrary case of diagonal disorder applies 
when the particles are positioned periodically but have inhomogeneous shapes or different 
resonance frequencies. For the sake of brevity, we limit our discussion to the case of pure 
diagonal disorder. All the techniques used in this article are also applicable in presence of 
off- diagonal disorder. 




(7) 




(8) 
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In the units described before, the off-diagonal elements of A4 are written as 



D. 



k\i-i\ {k\%-j\Y' 



1 I 1 



) 




(9) 



2 



(k\t-j\y {k\t-j\yj 




(10) 



for TEM and TM excitations, respectively. 

Since the Ohmic damping is absent and the lowest order radiation damping is independent 
of the particle geometry, the diagonal elements are inhomogeneous only in their real parts. 
We choose the real part from the the set of random numbers U{—W/2, W/2) which has a box 
probability distribution around zero with a width W. The imaginary part of the diagonal 
elements is constant in these units and equals —2t/3. Considering the linear dependence 
of the inverse of polarizability ([8]) on the particle volume and detuning from resonance 
frequency, realizing a uniform distribution is practical. 

D. Hypothetic models 

The results of the perturbation approximation disagrees with some of the trends observed 
in our numerical results for weakly coupled systems. To shed light on the origin of these 
observations, we have performed similar statistical analysis on extra hypothetical models. 
In these four models, step by step, we transform our model for TEM excitation to an 
ensemble of orthogonal random matrices, for which extensive results have been reported in 
the literature (see Ref. |3|] for a recent review). In all these ensembles the diagonal elements 
are real random numbers selected from the set U{— W/2, W/2). The distinction is in the 
off-diagonal elements which are defined as follows: 

HO, The matrices in this model are orthogonal and they are the closest to the frequently 
used PLRBM ensemble. The offdiagonal elements are random real numbers given by 



where hij is a randomly chosen from U{—1, 1); i.e. uniformly distributed in [-1,1]. 
HI, These matrices are the Hermitian counterpart of the TEM coupling matrix with a 
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randomized phase factor for each element: 



Df<', = A^e^"^-, (12) 

nHi _ f-)± -t<f>ji 
^i>j — ^ij^ ) 

where 0j<j is a random number from U{—7t, it). 

CI, This ensemble of complex-symmetric matrices resemble the TEM model with a ran- 
domized coupling phase. 

Dg.i ^ D^^e^'^^^, (13) 
where (pij = (pji are random numbers from f/(— vr, it). 

H2, This model is based on the Hermitian form of TEM interaction and the phase factor 
is kept periodically varying. 

Dfi ^ A^, (14) 

^i>j — ^ly 



III. ANALYTIC RESULTS 

Decomposition ([6]) relates the overall statistical behavior of the system to the properties 
of the eigenmodes and their corresponding eigenvalues. The dipole chain is an open system 
and the excitations are subject to radiation losses, which lead to the exponential decay 
of a mode. Therefor it is not possible to distinguish between disorder and loss origins of 
localization only based on the spatial extent of a mode. For these types of systems, statistical 
analysis has shown to be the only unambiguous method of studying Anderson transition. 
Therefor we study the scaling behavior. For this analysis, based on the eigenvectors in 
the position basis, two important indicators are considered: 1-the probability distribution 
function (PDF) of the wavefunction intensities and 2-the generalized inverse participation 
ratios (GIPR). 

The PDF is more easily accessible in experiments^. For numerical analysis, it has proven 
to be an accurate tool for measuring the scaling exponent in a finite size system^*^ and 
extracting the critical exponent from finite size scaling analysis^i. With the parametrization 
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V{a\ W, L, b) the PDF is sufficient for characterizing an Anderson transition. Here, a = 
In Ib/ ln{b/ L) with Ib = Yl^i=i l^n(3^j)l^ the integrated intensity over any box selection of 
length b. The effective disorder strength is parameterized by W, but the exact definition 
depends on the model. Criticality of eigenfunctions demands the scale invariance of the 
PDF. It means that the functional form of V does not change with system size for a fixed 
b/L. Away from the transition point, the maximum of the PDF, am, exhibits finite size 
scaling behavior—. This maximum shifts to higher(lower) values at the localized(extended) 
side of the transition. 

Another widely used set of quantities for evaluating the scaling exponents is the set of 
GIPR, which are proportional to the moments of the PDF. For each wavefunction GIPR are 
defined as 

L 
i=l 

At criticality, the ensemble averaged GIPR, (Pq), scales anomalously with the length L as 

(P,)~L-^^(^-^), (16) 

where dq is called the anomalous dimension. For multifractal wavefunctions, which are 
characteristic of Anderson transitions dq is a continuous function of q. From the definition. 
Pi = 1 and Pq = L. In practice, the GIPR can also be evaluated by box-scaling for a single 
system size, given a large enough sample^. 

A. Perturbation results for the weak coupling regime 

In the regime of weak coupling Wk ^ 1 the off-diagonal matrix elements of the Hamilto- 
nian are small compared to the diagonal ones. Therefore the moments of the eigenfunctions 
can be computed perturbatively using the method of the virial expansionr^'2^"— . To this 
end we generalize the route suggested in Ref. |2(|] to the case of the non-Hermitian random 
matrices. 

By using this perturbation analysis, we find that TEM eigenfunctions scale critically with 
the length of the system. The criticality is set by the inverse linear interaction term in Eq. ([9]) 
which dominates at large distances. In the weak-coupling regime, the set of multifractal 
exponents can be explicitly calculated. The result is different from the universal one found 
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for all critical models with Hermitian random matrices*^ and is given by 

2co(g) 1 

The detailed derivation of this result as well as an explicit expression for co{q) are pre- 
sented in Appendix |X1 The corresponding result for the orthogonal matrices reads 

_ 40Fr(g-l/2) 1 
- WkViq) ' ^ > 2- 

If a similar analysis is performed on the TM eigenfunction, the GIPR converge at large 
system sizes implying that the eigenf unctions are localized. This is due to the behavior 
of the coupling at large distances. 

In the following section, an extensive comparison is made between these analytical ex- 
pressions and the numerical simulations. 



IV. NUMERICAL RESULTS 



By direct diagonalization of a large ensemble of matrices, we have studied the PDF and 
GIPR scaling of the eigenfunctions of matrices from all the models introduced in the previous 
sections. Several values of disorder strength W and carrier wavenumber k are considered 
for matrices with sizes from L = 2^ to 2^^. Each matrix is numerically diagonalized with 
MATLAB using the ZGGEV algorithm. The number of analyzed eigenfunctions for each 
set of parameters is around 10^. Computation time for diagonalization of the largest matrix 
is 20 minutes on a PC. 



A. Spectrum of the homogeneous chain 

We start by analyzing the spectrum of the homogenous chain on resonance (W=0) where 
all the diagonal element are given hj Aia = — 2t/3. Typical spectra for TEM and TM 
excitations are shown in Figures [U^a) and[2t^a) for A; = 1. 

For k < 1.4, the TEM eigenvalues are divided into almost-real and complex subsets. The 
almost-real (Ime <^ Ree) subset corresponds to subradiative (bound) eigenstates. These 
eigenstates have a wavelength shorter than the free space propagation^!^ and cannot couple 
to the outgoing radiation, except at the two ends of the chain. The eigenmodes corresponding 
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to complex eigenvalues (line ~ Ree) are superradiative. For these modes a constructive 
interference in the far-field enhances the scattering from each particle in comparison with 
the an isolated one. 

From the form of expansion it is clear the eigenstates with (close to) zero eigenvalues 
will dominate the response of the system to external excitation. However, different regions 
in the spectrum can be experimentally probed by two approaches: Firstly, by changing the 
lattice spacing, or secondly, by detuning from the resonance frequency, which will add a 
constant real number to the diagonal of the interaction matrix Ai. Close to the resonance 
this number is linearly proportional to the frequency variation. This shift results in driving 
a different collective excitation, which has obtained the closest eigenvalue to the origin of 
the complex plane. 

B. The effect of disorder 

As mentioned before, disorder is introduced to the system by adding random numbers 
from the interval [—W/2, W/2] to the diagonal of A4. With this setting, the parameter space 
has two coordinates W and k, and g = {Wk)~^ is the coupling parameter. The weak and 
strong coupling regimes correspond to g <^ 1 and g ~ respectively. For A; < 0.5 the 
short range behavior is dominated by the quasi-static part of the interaction. The eigenstates 
of the disordered chain are thus exponentially decaying -similar to localized states. Since 
we are mainly interested in the critical behavior of eigenfunctions we focus on the region 
with 0.5 < A; < 3. 

1. Small disorder 

In the intermediate and strong coupling regime, the spectrum of the disordered matrix 
keeps the overall form of the homogeneous case (where the real part of the diagonal is zero), 
as can be seen in Figures [U^b) and|2]^b). 

At the sub-radiative band-edge of the TEM spectrum, modes of different nature mix due 
to disorder. This region is magnified in the inset of Fig.[T]^b). The eigenmodes corresponding 
to this region are of hybrid character. They consist of separate localization centers that are 
coupled via extended tails of considerable weight. The typical size of each localized section 
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is longer than the interparticle spacing. Similar modes have been observed in a quasi-static 
investigation of two-dimensional planar composites^I. For one-dimensional systems they are 
sometimes called necklace states in the literature^i^. Heuristically, this behavior can be 
attributed to the disorder induced mixing of sub-radiative and super-radiative modes which 
have closeby eigenvalues in the complex plain. Further evidence for this mixed behavior will 
be later discussed based on the shape of PDF in section IIV C 2[ 

For TM polarization all of the eigenstates become exponentially localized with power- 
law decaying tails. The localization length increases towards the band center. Therefore, 
in a chain with finite length, one will see two crossovers in the first Brilluoin zone, from 
localized to extended and back. However, the nature of localization seems to be different 
at the two ends. The subradiative modes [Ime <^ Ree) are localizaed due to interference 
effects similar to the Anderson localization while the superradiative modes (Im e ~ Re e) are 
localized by radiation losses. These two crossover regions eventually approach each other 
and disappear as the amount of disorder is increased, leading to a fully localized spectrum of 
eigenmodes. The spectral behavior is more complicated for higher wavenumbers with k > it 
but a discussion on that is further than the scope of this article. 

2. Large disorder 

In the weak coupling regime, the matrix is almost diagonal and thus the eigenvalues just 
follow the distribution of the diagonal elements. Typical eigenstates are shown in Fig. |3l 
As will be shown later, for this regime, all the eigenstates for TEM and TM are localized 
(since the coupling is weak) except for a band (about 20% width) of TEM eigenstates with 
the most negative real part of their eigenvalues. The states in this band show multifractal 
(critically extended) behavior for any arbitrarily weak coupling. Existence of these states is 
one of the major results of this investigation and their statistical analysis is the main subject 
of interest in the rest of this report. 

The multifractality of eigenfunction in the weak coupling regime is inline with the pre- 
diction of the virial expansion result f lT7|) . However our theory cannot describe why only 
a part of the TEM eigenstates are critically extended and the rest of them are localized, 
according to the numerical results. 
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C. Scaling behavior of PDF 



The scaling of PDF is an effective tool for analyzing the localized to extended transition 
in sample with finite Iength22i2ii22. We also use this statistical indicator to distinguish the 
regions of critical scaling. Only those eigenmodes for which their scaled PDF for different 
system sizes overlap are critical. For the wavefunctions that fulfill this criteria, the scaling 
of GIPR is analyzed. This second analysis confirms the presence of critical behavior by 
checking the the power-law scaling behavior. The logarithmic slope gives the multifractal 
dimensions. We have preformed extensive survey of the size-scaling behavior of PDF over 
the W, k space with 10^ eigenfunctions for each configuration. 

For each system size the scaled PDF is approximated by a histogram P(ln/5/ln(6/L)) 
over the sampled eigenfunctions These histograms are shown in Figures H] to [7] for different 
models. The shift of the peak of the distribution toward larger values (higher density of 
darker points) by an increase in the system size is a signature of eigenmode localization. A 
shift in the opposite direction towards a Guassian distribution with a peak at = d = 1 
is characteristic of the extended states. Overlap of these histograms signifies the critical 
behavior of the eigenmodes. 

1. TM and TEM modes in the weak coupling regime 

The typical scaling of PDF for TM modes is shown in Fig. |H It clearly reveals the 
localized behavior of these eigenfunctions. This is the generic behavior observed for these 
modes at any point in the parameter space. This result is in agreement with the Levitov's 
prediction, since the coupling is decaying as r~^. Localization in disordered one dimensional 
systems has already been studied extensively and we do not discuss it further in this article. 

In the regime of weak coupling, Wk > 10, the numerical results show convincing indica- 
tion of critical scaling in a band of the TEM modes. These results are plotted in Fig. [5l The 
band of critical modes consists of those with the most negative real part of their eigenvalues. 
Outside this band the eigenfunctions show scaling behavior similar to localized modes as 
shown in Fig. [5]^a) and (b). This crossover from localized to critical eigenfunctions may be 
useful for measuring the critical exponent. However, critical exponent must be defined based 
on a proper ordering of eigenvalues, which is known to be a non-trivial task for complex 
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eigenvalues. 



2. TEM modes in the strong coupling regime 

In the strong coupling regime, i.e. weak disorder, we have found it more representative 
to order the complex eigenvalues by their argument. A narrow region near the negative real 
axis is selected as shown in Fig. [UJ^b). The histograms representing the scaled PDF of these 
eigenf unctions are plotted in Fig. El^a). These histograms do not overlap so the criticality 
cannot be verified. Meanwhile, the behavior is neither representative of the localized modes 
nor the extended modes. It appears that the overall extent of the state is comparable 
with the system size even for the longest chain, but it has an strongly fluctuating internal 
structure, similar to critical states. Typical eigenmodes of this regime are shown in Fig.[6](b). 
Since the scaled PDF histograms do not overlap, we cannot prove the multifractal nature of 
the states with a formal logarithmic scaling. Describing the true nature of these modes and 
their statistical behavior needs further theoretical modeling. 

3. PDF of the intermediate models 

The results of the perturbation calculations in Sec. IIII Al are insensitive to the details of 
the model. Therefore they cannot describe some of our observations that are based on direct 
numerical diagonalization. For example, according to the perturbation theory, all the TEM 
modes in the weak coupling regime must be critical. This prediction does not agree with 
the simulation results since PDF scaling in observed for only part of these modes. The same 
numerical analysis on Hermitian random banded matrices perfectly matches the results of 
perturbation results. 

To further explore the origin of this deviation for complex-symmetric matrices of our 
model for TEM excitations, we have performed the same numerical procedures on the hypo- 
thetic models introduced in Sec. IIIDi The PDF scaling graphs for these models are depicted 
in Fig. [71 All these results are for the regime of weak coupling with the same W and k. 

HO, The matrices in this model are orthogonal and they are the closest to the frequently 
used PLRBM ensemble with an interaction decay exponent /i = 1. The PDF shows 
perfect scaling as depicted in Fig.[71^a). The statistics is obtained by sampling from 12% 
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of the eigenvectors at the band center, with eigenvalues closest to zero. The analysis 
shows the same critical behavior (not shown) for the two ends of the spectrum. These 
results also confirm that our choice of numerical precision and sampling is sufficient 
for the essential conclusions we get. 

HI, These matrices are also Hermitian like model HO. The magnitude of the off-diagonal 
elements is not random, but follows the decay profile of TEM complex-valued cou- 
pling on]). Only the phase is randomized. Critical scaling of the eigenf unctions is again 
evident from the PDF scaling depicted in Fig. UIJo)- 

CI, This ensemble of complex-symmetric matrices resembles the TEM model. The phases 
of the off-diagonal elements are randomized like the model HI. The finite-size scaling 
of PDF, depicted in Fig. Wic) shows the behavior that is attributed to localized modes. 
For localized eigenvectors the peak of the distribution shifts toward higher values of 
a, which signifies a higher density for points with a low intensity. 

H2, This model is the Hermitian form of TEM coupling matrix. The difference between 
this model and HI is in the phase factor, which is kept periodic like the original Green 
function. The only random elements of these matrices are the diagonal ones. Despite 
the minor difference between models HI and H2, the result of PDF scaling analysis is 
completely different. These results are depicted in Fig. Wl^d) and show that the eigen- 
vectors are localized. This observation is inconsistent with the perturbation theory, 
which predicts critical behavior for this model like HO and HI. It is important to point 
out that the considered periodicity for the interaction phase /c = 1 is incommensurate 
with the periodicity of the lattice, which equals 27r in our redefinition of units. 

D. Multifractal analysis 

Since the critical scaling of part of the TEM eigenmodes in the weak coupling regime 
is clearly observed in the scaling of PDF, we apply generic techniques of multifractal (MF) 
analysis to quantify the MF spectrum and compare it with our theoretical results. 

We have used both size scaling and box scaling methods for extracting the scaling ex- 
ponents of GIPR for several different parameters. We do not observe significant differences 
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in the results of either method (comparison not shown). Therefore, due to its faster com- 
putation, we use the box scahng analysis on the largest system sizes, L = 4096, to extract 
the anomalous exponents dg for several values of W and k. A summary of these results for 
different values of disorder strength is depicted in Figures [S] and M To show the precision of 
the numerical analysis, we have also performed this analysis for Hermitian model HI. The 
results are shown in Fig. [8](b) and compared with the theoretical prediction of Eq. (ITSl) . Ex- 
cellent matching between theory and simulation is evident for the Hermitian case. However, 
for the complex-symmetric matrices (corresponding to TEM coupling) the numerical results 
show significant deviations from the prediction of perturbation analysis, indicating that the 
first order virial expansion is insufficient for describing that model. 

In particular, according to Eq. f[T7|) dg is proportional to the coupling strength g = (Wk)~^ 
in the weak coupling regime. The results of direct diagonalization show, in contrast, a 
dependence of dg on 14^ at a fixed value of g. This fact can be seen in Fig. [8]^a). The 
numerical results are systematically lower than the theoretical prediction for k < 3. 

Furthermore, the dependence of MF dimensions on the coupling strength is checked for 
9 < Wk < 150. The results are shown in Fig. |9]for k = 1 and k = 3. The overall inverse 
linear behavior is observed for Wk > 30. But the quantitative correspondence between the 
numerical results and prediction ( [T7|) from perturbation analysis is only met for the large 
values of k and high moments of GIPR, g > 3. 

E. The singularity spectrum 

Another frequently used representation of the multifractality is called the singularity 
spectrum, f{a). This representation is completely analogous to the one using anomalous 
dimensions. For the sake of completeness, we also show this representation for two of the 
models that are critical in their scaling behavior. 

The singularity spectrum f{a) is the fractal dimension of the subset of those points in 
the wavefunction for which intensity scales as It is related to the set of anomalous 

exponents r(g) = dg{q — 1) by a Legendre transform 

f(a) =qa-T{q), a = T'{q) q = f'{a),. (19) 

The quantity a introduced here is related by an irrelevant scaling prefactor to a, which was 
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used for the definition of PDF. However, specially for a skewed PDF, this prefactor can 
significantly deviate from unity and therefore a and a are not exactly equivalent quantities. 

For a precise derivation of a and /(a) one has to perform a full scaling analysis on the 
intensity distribution. This is either possible by applying relation f lT^ to the calculated 
set of anomalous exponents or by a direct processing of the wavefunction intensities. The 
latter method, which was introduced by Chhabra and Jensen^i, is computationally superior. 
For this method there is no need for a Legendre transform, which is very sensitive to the 
numerical uncertainties. 

We have applied the direct determination method to extract the singularity spectrum for 
TEM critical eigenfunctions and the Hermitian model HO. The results are shown in Fig. (TUi 
As can be seen in both graphs, the position of the peak of the spectrum is different from 
the peak of the corresponding PDF plots, (5^, which are presented in Figures |5t^c) and[7l^a). 
This difference is due to the large skewness of the PDF resulted from the very weak coupling 
regimes that are considered in this article. 

For the Hermitian critical models the domain of a is restricted to (0, 1(£) due to the 
symmetry relation^ of /(a). Yet, there is no proof that this symmetry also holds for non- 
Hermitian matrices. From our data, it seems plausible that this symmetry is actually broken 
and there are some points with a > 2. However, to provide a strong numerical evidence for 
this statement one has to analyze much larger ensembles with higher numerical precision, 
which is beyond the scope of the current article. 

V. SUMMARY AND CONCLUSION 

We have investigated, theoretically and numerically, the statistical properties of the eigen- 
modes of a class of complex-symmetric random matrices, which describe the electromagnetic 
propagation of polarization waves in a chain of resonant scatterers. We have found that all 
of the TM modes are localized in the weak coupling regime. The TEM modes in this regime 
show critical behavior due to the dependence in the dyadic Green function. This crit- 
ical behavior is in agreement with the results of the method of virial expansion for almost 
diagonal matrices. We have used this method to calculate the MF spectrum of TEM modes. 

Although the perturbation theory suggests criticality for all TEM modes, the numerical 
analysis shows this type of scaling only for part of the spectrum in the complex plain. This 
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is understandable in the sense that the first oder result of the perturbative approach gives an 
oversimplified picture, which is insensitive to details of the model such as a non-trivial phase 
dependence of the matrix elements. To reveal which aspect of the TEM coupling accounts for 
the existence of a critical band in the spectrum, we have analyzed three intermediate models. 
These models have properties between the dipole chain interaction matrix and power-law 
Hermitian banded random matrices. The summary of the scaling results for all these models 
is shown in Fig. [11] It seems that both non-Hermitian character of the TEM coupling and 
the periodic phase of the interaction between dipoles is important for the observed critical 
eigenmodes. 

Our analysis also resulted in another unexpected finding. The eigenvectors of Hermitian 
banded matrices with coupling are no longer critically scaling if the interaction phase 
is set periodically. In our model H2, the randomness is only on the diagonal. Based on the 
PDF scaling results, we clearly see that the eigenvectors are localized. This is in contrast 
with the commonly believed conjecture that an interaction potential with a phase that is 
incommensurate with the lattice can be considered as random. 

Criticality of wavefunctions has been studied theoretically and numerically for several 
models in the context of condensed matter physics. Recently, such wavefunctions have been 
observed near the Anderson transition for elastic waves^^ and electronic density of states 
at an interface^. The recent advances in optical and microwave instrumentation makes 
it possible to experiment in details the propagation of electromagnetic waves in artificially 
made structures. Our report points out to those systems in which such critical phenomena 
can be directly measured. These measurements provide a lot of insight for generic models 
of wave transport in disordered system. 
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Appendix A: Perturbation results for eigenvector moments 



In the regime of the strong multifractahty, i.e. when Wk ^ 1, the density of states is 
determined only by the distribution of the diagonal elements and hence is uniform. Therefore 
we assume that all eigenvectors are characterized by the same scaling exponents and we 
define the ensemble average GIPR as 

n=l i=l 

The moments of the eigenvectors can be extracted from the powers of the diagonal elements 
of the Green functions and therefore can be computed using the method of the the virial 
expansion. 

In the first order of the virial expansion, which corresponds to a pure diagonal matrix, the 
eigenvectors consists of only one non-zero component and therefore all moments are equal 
to one due to normalization: (-Pq)^^^ = 1. 

In the next order of the virial expansion the contribution to the eigenvectors moments 
from all possible pairs of two levels of the unperturbed system should be taken into account. 
Thus we need to calculate the moments of the eigenvectors of the following 2x2 matrices: 

M(U) = ( . (A2) 

Denoting by j) the moments of the second component of the corresponding eigenvectors 
we have 

{P,)^'^ = \Y.{PS,3)-1) (A3) 

The subtraction of 1 in the equation above eliminates the contribution already taken into 

account in the diagonal approximation. 

Let us introduce the following notation for M{i.j): 

. El he'' , 

= I . h (A4) 



2 



where Ei = Ei, E2 = Ej and /le*"^ = Mij — Mji- Below we also denote Pq{i,j) by Pg. As 
the eigenvectors of H{i,j) do not depend on the absolute values of the matrix elements, but 
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only on their relative ratios, we may assume now that Ei and E2 are distributed uniformly 
in [—1/2, 1/2] and h = e^''^^~^^ /Wk\i — j\- For the reasons which will be described later, the 
other terms in the coupling elements are not considered here. By the assumption of weak 
coupling \h\ ^ 1 and writing the eigenvectors of M{i,j) explicitly, we obtain: 



P. 



Q{x) 



f'^ dE, f'^ dE^iQUE, - E2)/h) + Q{{E2 - E,)/h)), 

1-1/2 J-1/2 

q 

4 



(A5) 



In the limit h ^ one can show that Q{r) = 9(r) and thus Pq = 1. This result corresponds 
to the diagonal approximation and is cancelled by —1 in Eq.( 1A3p . In order to find first 
non-trivial contribution we need to compute a term which is linear in h. To this end we first 
differentiate Eq. (]A5|) with respect to h, then change the integration variables {Ei,E2} — )■ 
{E\, X = {El — E2)/h] and consider the limit /i — )■ 0: 

dPr. 



1™ Jl. 

h^o dh 



The expansion of P„ then takes the form: 



dx xQ'{x) = —F{q, 



P, = l-\h\Fiq,<f^) + Oih'). 
Collecting the contributions from all the off-diagonal elements we obtain 



1 1 
WkL 



F{q,k\i-j\) 



Expanding F{q, 0) in the Fourier series F{q, 0) = J^p^^pil)^''^'^ 



1 



1 co(g) 



1 ^ ^ Cp(g)e*P'=l*-^l 



2co(g)lnL + 0(l), 



so that we arrive at the following result: 



2co(g) 



InL, co(g) 



1 



27r 



dx xQ'J^x, g. 



where Q is defined in Eq. OASp and can be written as 

4 



Q(a;,g, 



4 + 



X 



- \fxFV^e-^ 



(A6) 



(A7) 



(A8) 



(A9) 



(AlO) 



(All) 
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Comparing this result with the scahng of the moments by changing the system size 

(Pg) oc L-'^''(«-i), (A12) 



we find the expressions for the fractal dimensions 

One can check that in the case of the orthogonal matrices, i.e. when = 0, the old result 

can be reproduced after the change of the variable x = {2w — 1)/ \/w{l — w): 



/oo 
dx xQ'^{x,q,0) 
-oo 

^ 2 dw- ^ ' 



w 



The result for complex-symmetric matrices deviate only slightly from the one for orthogonal 
case. This deviation is most pronounced around q — 2.5. 
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FIG. 1: Complex valued spectrum of the TEM interaction matrix ([9]) with k = 1 for (a) homoge- 
neous {W = 0) and disordered in regimes of (b) strong {W = 2) and (c) weak {W = 20) coupling. 
The dashed square in (c) shows the region where the corresponding TEM eigenmodes are scaling 
critically. The eigenmodes corresponding to the eigenvalues in the dashed regions are selected for 
further statistical analysis. 
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FIG. 2: Same as Fig. [T]for TM polarization. 
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FIG. 3: Typical TEM critical (red solid line) and TM localized (black dashed line) eigenvectors 
of matrices defined in section III CI with k = 1 and W = 10. The inset shows the same plot in 
logarithmic scale. 




FIG. 4: Scaling of PDF for TM eigenmodes for different lengths of the chain and correspondingly 
scaled box sizes of 6 = x L. Different line types (and colors) correspond to different system 
sizes as indicated by the legend. = 30 and k = 1. The shift of the peak to the larger values of 
a indicates that the eigenmodes are localized. 
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FIG. 5: Same as Fig.[4]for TEM eigenmodes. For each figure 12% of the eigenmodes are used with 
(a) most positive, (b) closest to zero, and (c) most negative real part of their eigenvalues. Critical 
scaling is only observed in (c). The shift in the peak of the distribution shows that the rest of the 
eigenmodes are localized. 
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FIG. 6: (a) Similar to Fig. [H for TEM modes in the regime of strong coupling with k = 1 and 
W = 0.7. The analyzed eigenmodes are selected from the spectral region indicated by the dashed 
triangle in Fig. (D^b). (b) Typical eigenmodes used for the PDF in (a) for two system lengths 
L = 4096 and L = 512 (inset). 
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FIG. 7: Similar to Fig. |4]for models (a) HO, (b) HI, (c) H2, and (d) CI. Critical scaling is observed 
for HO and HI models. The eigenvectors of the H2 and CI models are localized. These models are 
defined in Sec. IHDI 
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FIG. 8: The anomalous dimensions dq (symbols) at a fixed coupling strength for (a) TEM critical 
modes and (b) Hermitian matrices HI are compared with the corresponding results (solid line) 
(|17|) an d (1181) from perturbation analysis. For Hermitian matrices, the symmetry relation predicted 
in Ref. 84] is used for plotting the theoretical curve at negative q. The errors estimated from the 
least squares fitting routine are smaller than the symbol sizes and are not shown. The largest error 
in dg for point q = 5.5 on the graph is ±0.02. 




FIG. 9: The anomalous dimensions dg extracted from direct numerical diagonalization (symbols) 
for three different values of q are compared with the perturbation results (lines) of Eq. (I17p in the 
weak coupling regime for different values of disorder and (a) A; = 1 or (b) A; = 3. The numerical 
results converge to the theory slowly. 
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FIG. 10: The multifractal spectrum f{a) for (a) TEM critical modes and (b) Hermitian matrices 
HO extracted directly from the eigenvectors by using the method of Chhabra and Jensen™. For 
both graphs W = 30 and k = 1. The error bars indicate to the standard deviation among 20 
reahzations of disorder and are smaher than the symbol size for most of the data points. The 
dashed lines are guides to the eye. 
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FIG. 11: Summary of the scaling analysis on PDF of eigenvectors of matrices from various models. 
The intermediate hypothetic models transform the Hermitian REM to the model describing TEM 
coupling. The colored boxes indicate those models that show critical scaling behavior. The model 
indicated by white boxes have localized eigenvectors. For TEM modes, only a part of the spectrum 
is critical. 
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